# total random sites to create
tot <- nrow(restdat)
id <- stri_rand_strings(tot, 4)
dts <- range(restdat$date)
# rnorm(10 * tot, mean(reststat$lon), sd(reststat$lon))
# rnorm(10 * tot, mean(reststat$lat), sd(reststat$lat))
ext <- bbox(tbpoly)
lon <- runif(10 * tot, ext[1, 1], ext[1, 2])
lat <- runif(10 * tot, ext[2, 1], ext[2, 2])
tmp <- SpatialPoints(cbind(lon, lat),
proj4string = crs(tbpoly)
) %>%
.[tbpoly, ] %>%
.[sample(1:nrow(.@coords), tot, replace = F), ]
restdat_rnd <- tibble(id) %>%
mutate(
date = sample(seq(dts[1], dts[2]), tot, replace = T),
top = sample(c('hab', 'wtr'), tot, replace = T)
)
reststat_rnd <- tibble(id) %>%
mutate(
lat = tmp$lat,
lon = tmp$lon
)
resgrp <- 'top'
restall_rnd <- left_join(restdat_rnd, reststat_rnd, by = 'id')
names(restall_rnd)[names(restall_rnd) %in% resgrp] <- 'Restoration\ngroup'
# base map
ext <- make_bbox(reststat_rnd$lon, reststat_rnd$lat, f = 0.1)
map <- get_stamenmap(ext, zoom = 10, maptype = "toner-lite")
pbase <- ggmap(map) +
theme_bw() +
theme(
axis.title.x = element_blank(),
axis.title.y = element_blank()
)
# map by restoration type
pbase +
geom_point(data = restall_rnd, aes(x = lon, y = lat, fill = `Restoration\ngroup`), size = 4, pch = 21)
# map by date
pbase +
geom_point(data = restall_rnd, aes(x = lon, y = lat, fill = factor(date)), size = 4, pch = 21)
# barplot of date counts
toplo <- restall_rnd %>%
group_by(date)
ggplot(restall_rnd, aes(x = factor(date))) +
geom_bar() +
coord_flip() +
theme_bw() +
theme(
axis.title.y = element_blank()
) +
scale_y_discrete(expand = c(0, 0))
wqmtch_rnd <- get_clo(restdat_rnd, reststat_rnd, wqstat, resgrp = 'top', mtch = mtch)
save(wqmtch_rnd, file = 'data/wqmtch_rnd.RData', compress = 'xz')
head(wqmtch_rnd)
## # A tibble: 6 x 5
## stat resgrp rnk id dist
## <int> <chr> <int> <chr> <dbl>
## 1 47 hab 1 06R7 3863.373
## 2 47 hab 2 8o2A 4035.380
## 3 47 hab 3 b5PN 5989.463
## 4 47 hab 4 NRgK 7774.976
## 5 47 hab 5 W146 9705.679
## 6 47 hab 6 gTaE 9821.184
##
# plots
# combine lat/lon for the plot
toplo <- wqmtch_rnd %>%
left_join(wqstat, by = 'stat') %>%
left_join(reststat_rnd, by = 'id') %>%
rename(
`Restoration\ngroup` = resgrp,
`Distance (dd)` = dist
)
# restoration project grouping column
resgrp <- 'top'
restall_rnd <- left_join(restdat_rnd, reststat_rnd, by = 'id')
names(restall_rnd)[names(restall_rnd) %in% resgrp] <- 'Restoration\ngroup'
# extent
ext <- make_bbox(wqstat$lon, wqstat$lat, f = 0.1)
map <- get_stamenmap(ext, zoom = 12, maptype = "toner-lite")
# base map
pbase <- ggmap(map) +
theme_bw() +
theme(
axis.title.x = element_blank(),
axis.title.y = element_blank()
) +
geom_point(data = restall_rnd, aes(x = lon, y = lat, fill = `Restoration\ngroup`), size = 4, pch = 21) +
geom_point(data = wqstat, aes(x = lon, y = lat), size = 2)
# closest
toplo1 <- filter(toplo, rnk %in% 1)
pbase +
geom_segment(data = toplo1, aes(x = lon.x, y = lat.x, xend = lon.y, yend = lat.y, alpha = -`Distance (dd)`, linetype = `Restoration\ngroup`), size = 1)
# closest five percent
fvper <- max(toplo$rnk) %>%
`*`(0.2) %>%
ceiling
toplo2 <- filter(toplo, rnk %in% c(1:fvper))
pbase +
geom_segment(data = toplo2, aes(x = lon.x, y = lat.x, xend = lon.y, yend = lat.y, alpha = -`Distance (dd)`, linetype = `Restoration\ngroup`), size = 1)
# closest all combo
toplo3 <- toplo
pbase +
geom_segment(data = toplo3, aes(x = lon.x, y = lat.x, xend = lon.y, yend = lat.y, alpha = -`Distance (dd)`, linetype = `Restoration\ngroup`), size = 1)
Get weighted average of project type, treatment (before, after) of salinity for all wq station, restoration site combinations.
salchgout_rnd <- get_chg(wqdat, wqmtch_rnd, statdat, restdat_rnd, wqvar = 'sal', yrdf = yrdf, chgout = TRUE)
salchg_rnd <- get_chg(wqdat, wqmtch_rnd, statdat, restdat_rnd, wqvar = 'sal', yrdf = yrdf)
save(salchgout_rnd, file = 'data/salchgout_rnd.RData')
save(salchg_rnd, file = 'data/salchg_rnd.RData')
head(salchgout_rnd)
## # A tibble: 6 x 3
## # Groups: stat [2]
## stat cmb cval
## <int> <chr> <dbl>
## 1 6 hab_aft 23.62066
## 2 6 hab_bef 23.45522
## 3 6 wtr_aft 24.24440
## 4 6 wtr_bef 24.26601
## 5 7 hab_aft 23.80824
## 6 7 hab_bef 24.33082
head(salchg_rnd)
## # A tibble: 6 x 4
## stat hab wtr cval
## <int> <fctr> <fctr> <dbl>
## 1 6 hab_aft wtr_aft 23.93253
## 2 6 hab_aft wtr_bef 23.94334
## 3 6 hab_bef wtr_aft 23.84981
## 4 6 hab_bef wtr_bef 23.86062
## 5 7 hab_aft wtr_aft 24.39695
## 6 7 hab_aft wtr_bef 24.11447
Get conditional probability distributions for the restoration type, treatment effects, salinity as first child node in network.
wqcdt_rnd <- get_cdt(salchg_rnd, 'hab', 'wtr')
head(wqcdt_rnd)
## # A tibble: 4 x 5
## hab wtr data crv prd
## <fctr> <fctr> <list> <list> <list>
## 1 hab_aft wtr_aft <tibble [45 x 2]> <dbl [2]> <data.frame [100 x 3]>
## 2 hab_aft wtr_bef <tibble [45 x 2]> <dbl [2]> <data.frame [100 x 3]>
## 3 hab_bef wtr_aft <tibble [45 x 2]> <dbl [2]> <data.frame [100 x 3]>
## 4 hab_bef wtr_bef <tibble [45 x 2]> <dbl [2]> <data.frame [100 x 3]>
Discretization of salinity conditional probability distributions:
salbrk_rnd <- get_brk(wqcdt_rnd, qts = c(0.33, 0.66), 'hab', 'wtr')
salbrk_rnd
## # A tibble: 8 x 5
## hab wtr qts brk clev
## <fctr> <fctr> <dbl> <dbl> <dbl>
## 1 hab_aft wtr_aft 25.89943 0.4078706 1
## 2 hab_aft wtr_aft 29.44066 0.8618699 2
## 3 hab_aft wtr_bef 25.59620 0.3693579 1
## 4 hab_aft wtr_bef 29.38185 0.8551371 2
## 5 hab_bef wtr_aft 26.36323 0.4313804 1
## 6 hab_bef wtr_aft 29.94623 0.8684977 2
## 7 hab_bef wtr_bef 26.24658 0.4148461 1
## 8 hab_bef wtr_bef 29.98210 0.8702296 2
A plot showing the breaks:
toplo <- dplyr::select(wqcdt_rnd, -data, -crv) %>%
unnest
ggplot(toplo, aes(x = cval, y = cumest)) +
geom_line() +
geom_segment(data = salbrk_rnd, aes(x = qts, y = 0, xend = qts, yend = brk)) +
geom_segment(data = salbrk_rnd, aes(x = min(toplo$cval), y = brk, xend = qts, yend = brk)) +
facet_grid(hab ~ wtr) +
theme_bw()
Get conditional probability distributions for the restoration type, treatment effects, salinity levels, chlorophyll as second child node in network.
# get chlorophyll changes
chlchgout_rnd <- get_chg(wqdat, wqmtch_rnd, statdat, restdat_rnd, wqvar = 'chla', yrdf = yrdf, chgout = TRUE)
chlchg_rnd <- get_chg(wqdat, wqmtch_rnd, statdat, restdat_rnd, wqvar = 'chla', yrdf = yrdf)
save(chlchgout_rnd, file = 'data/chlchgout_rnd.RData')
save(chlchg_rnd, file = 'data/chlchg_rnd.RData')
# merge with salinity, bet salinity levels
salbrk_rnd <- salbrk_rnd %>%
group_by(hab, wtr) %>%
nest(.key = 'levs')
allchg_rnd <- full_join(chlchg_rnd, salchg_rnd, by = c('hab', 'wtr', 'stat')) %>%
rename(
salev = cval.y,
cval = cval.x
) %>%
group_by(hab, wtr) %>%
nest %>%
left_join(salbrk_rnd, by = c('hab', 'wtr')) %>%
mutate(
sallev = pmap(list(data, levs), function(data, levs){
# browser()
out <- data %>%
mutate(
saval = salev,
salev = cut(salev, breaks = c(-Inf, levs$qts, Inf), labels = c('lo', 'md', 'hi')),
salev = as.character(salev)
)
return(out)
})
) %>%
dplyr::select(-data, -levs) %>%
unnest
salchg_rnd <- dplyr::select(allchg_rnd, stat, hab, wtr, salev, saval)
save(salchg_rnd, file = 'data/salchg_rnd.RData', compress = 'xz')
chlcdt_rnd <- get_cdt(allchg_rnd, 'hab', 'wtr', 'salev')
save(chlcdt_rnd, file = 'data/chlcdt_rnd.RData', compress = 'xz')
chlbrk_rnd <- get_brk(chlcdt_rnd, c(0.33, 0.66), 'hab', 'wtr', 'salev')
chlbrk_rnd %>%
print(n = nrow(.))
## # A tibble: 24 x 6
## hab wtr salev qts brk clev
## <fctr> <fctr> <chr> <dbl> <dbl> <dbl>
## 1 hab_aft wtr_aft lo 12.599474 0.4362990 1
## 2 hab_aft wtr_aft lo 17.177905 0.8734528 2
## 3 hab_aft wtr_aft md 6.709566 0.2432915 1
## 4 hab_aft wtr_aft md 8.343154 0.7019913 2
## 5 hab_aft wtr_aft hi 4.504318 0.4580072 1
## 6 hab_aft wtr_aft hi 5.367021 0.8636019 2
## 7 hab_aft wtr_bef lo 12.123996 0.3361189 1
## 8 hab_aft wtr_bef lo 16.457887 0.7989985 2
## 9 hab_aft wtr_bef md 8.037491 0.3559850 1
## 10 hab_aft wtr_bef md 10.811197 0.8275772 2
## 11 hab_aft wtr_bef hi 3.968360 0.3206374 1
## 12 hab_aft wtr_bef hi 4.668884 0.7478114 2
## 13 hab_bef wtr_aft lo 11.818910 0.3924014 1
## 14 hab_bef wtr_aft lo 16.028895 0.8215811 2
## 15 hab_bef wtr_aft md 6.693036 0.4236565 1
## 16 hab_bef wtr_aft md 8.007187 0.8770172 2
## 17 hab_bef wtr_aft hi 4.780222 0.4801252 1
## 18 hab_bef wtr_aft hi 5.800910 0.8786539 2
## 19 hab_bef wtr_bef lo 11.457732 0.3152265 1
## 20 hab_bef wtr_bef lo 15.537479 0.7633651 2
## 21 hab_bef wtr_bef md 8.090232 0.5215649 1
## 22 hab_bef wtr_bef md 10.561438 0.9094215 2
## 23 hab_bef wtr_bef hi 4.184045 0.3385945 1
## 24 hab_bef wtr_bef hi 4.999934 0.7635401 2
Final combinations long format:
chlbar_rnd <- chlbrk_rnd %>%
group_by(hab, wtr, salev) %>%
nest %>%
mutate(
data = map(data, function(x){
brk <- x$brk
out <- data.frame(
lo = brk[1], md = brk[2] - brk[1], hi = 1 - brk[2]
)
return(out)
})
) %>%
unnest %>%
gather('chllev', 'chlval', lo:hi) %>%
mutate(
salev = factor(salev, levels = c('lo', 'md', 'hi')),
chllev = factor(chllev, levels = c('lo', 'md', 'hi'))
)
save(chlbar_rnd, file = 'data/chlbar_rnd.RData', compress = 'xz')
chlbar_rnd %>%
print(n = nrow(.))
## # A tibble: 36 x 5
## hab wtr salev chllev chlval
## <fctr> <fctr> <fctr> <fctr> <dbl>
## 1 hab_aft wtr_aft lo lo 0.4362990
## 2 hab_aft wtr_aft md lo 0.2432915
## 3 hab_aft wtr_aft hi lo 0.4580072
## 4 hab_aft wtr_bef lo lo 0.3361189
## 5 hab_aft wtr_bef md lo 0.3559850
## 6 hab_aft wtr_bef hi lo 0.3206374
## 7 hab_bef wtr_aft lo lo 0.3924014
## 8 hab_bef wtr_aft md lo 0.4236565
## 9 hab_bef wtr_aft hi lo 0.4801252
## 10 hab_bef wtr_bef lo lo 0.3152265
## 11 hab_bef wtr_bef md lo 0.5215649
## 12 hab_bef wtr_bef hi lo 0.3385945
## 13 hab_aft wtr_aft lo md 0.4371538
## 14 hab_aft wtr_aft md md 0.4586999
## 15 hab_aft wtr_aft hi md 0.4055948
## 16 hab_aft wtr_bef lo md 0.4628796
## 17 hab_aft wtr_bef md md 0.4715923
## 18 hab_aft wtr_bef hi md 0.4271740
## 19 hab_bef wtr_aft lo md 0.4291797
## 20 hab_bef wtr_aft md md 0.4533607
## 21 hab_bef wtr_aft hi md 0.3985287
## 22 hab_bef wtr_bef lo md 0.4481385
## 23 hab_bef wtr_bef md md 0.3878566
## 24 hab_bef wtr_bef hi md 0.4249455
## 25 hab_aft wtr_aft lo hi 0.1265472
## 26 hab_aft wtr_aft md hi 0.2980087
## 27 hab_aft wtr_aft hi hi 0.1363981
## 28 hab_aft wtr_bef lo hi 0.2010015
## 29 hab_aft wtr_bef md hi 0.1724228
## 30 hab_aft wtr_bef hi hi 0.2521886
## 31 hab_bef wtr_aft lo hi 0.1784189
## 32 hab_bef wtr_aft md hi 0.1229828
## 33 hab_bef wtr_aft hi hi 0.1213461
## 34 hab_bef wtr_bef lo hi 0.2366349
## 35 hab_bef wtr_bef md hi 0.0905785
## 36 hab_bef wtr_bef hi hi 0.2364599
Discretesize chlorophyll data, all stations:
# discretize all chl data by breaks
chlbrk_rnd <- chlbrk_rnd %>%
group_by(hab, wtr, salev) %>%
nest(.key = 'levs')
allchg_rnd <- allchg_rnd %>%
group_by(hab, wtr, salev) %>%
nest %>%
full_join(chlbrk_rnd, by = c('hab', 'wtr', 'salev')) %>%
mutate(
lev = pmap(list(data, levs), function(data, levs){
out <- data %>%
mutate(
lev = cut(cval, breaks = c(-Inf, levs$qts, Inf), labels = c('lo', 'md', 'hi')),
lev = as.character(lev)
)
return(out)
})
) %>%
dplyr::select(-data, -levs) %>%
unnest %>%
rename(
chlev = lev,
chval = cval
)
save(allchg_rnd, file = 'data/allchg_rnd.RData', compress = 'xz')
A bar plot of splits:
ggplot(chlbar_rnd, aes(x = chllev, y = chlval, group = salev, fill = salev)) +
geom_bar(stat = 'identity', position = 'dodge') +
facet_grid(hab ~ wtr) +
theme_bw()
A plot showing the breaks:
toplo <- dplyr::select(chlcdt_rnd, -data, -crv) %>%
unnest %>%
mutate(
salev = factor(salev, levels = c('lo', 'md', 'hi'))
)
chlbrk_rnd <- unnest(chlbrk_rnd)
ggplot(toplo, aes(x = cval, y = cumest, group = salev, colour = salev)) +
geom_line() +
geom_segment(data = chlbrk_rnd, aes(x = qts, y = 0, xend = qts, yend = brk)) +
geom_segment(data = chlbrk_rnd, aes(x = min(toplo$cval), y = brk, xend = qts, yend = brk)) +
facet_grid(hab ~ wtr, scales = 'free_x') +
theme_bw()